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Abstract 

A discrete random walk method on grids was proposed and used to solve 
the linearized Poisson-Boltzmann equation (LPBE) Jjj]. Here, we present 
a new and efficient grid-free random walk method. Based on a modified 



"Walk On Spheres" (WOS) algorithm Q for the LPBE, this Monte Carlo 
algorithm uses a survival probability distribution function for the random 
walker in a continuous and free diffusion region. The new simulation method 
is illustrated by computing four analytically solvable problems. In all cases, 
excellent agreement is observed. 
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Random walk methods have been used to solve a wide variety of parabolic and elliptic 
partial differential equations (PDEs) []l]-f| . Generally, there are two broad classes of random 
walk methods; one uses discrete random walks on grids and the other continuous random 
walks in free space @-|5[]. One of the widely used continuous random walk methods, the 
"Walk On Spheres" (WOS) method |2|,|6|-|8[] , uses the first-passage probability distribution 
on a sphere to facilitate large steps in random walks. (The first-passage probability, w(x; x ), 
is the probability of hitting the vicinity of x on the bounding surface for the first time when 
the random walker starts from x , a point inside the bounding surface.) This continuous 
random walk method needs to discretize neither space nor time, nor the diffusing trajectory, 
and so it is particularly advantageous when the geometry of the region of interest is very 
complex or if the solution of the PDE is required at only a relatively small number of points. 

We are interested in solutions to the Dirichlet problem for the linearized Poisson- 
Boltzmann equation (LPBE) in the domain Q: 

vV(x) = «V(x), x g n, (i) 



V>(x) = Vo(x), x g an, (2) 

where k is called the inverse Debye length fl[]. Notice that when k 2 is zero, the above 
problem becomes a Dirichlet problem for the Laplace equation. The WOS method for the 
Dirichlet problem for the Laplace equation has been widely used P,|6|,|7|j3|j9|,p!0[1 . We will 
combine this WOS method with a survival probability density function which incorporates 
the term involving k 2 in the LPBE. 

In a discrete random walk method [Q] for the LPBE, the corresponding Master equation 
relates k 2 to the removal probability of the random walker on the grid. During each step 
of the discrete random walk, the walker either moves to one of the neighboring sites, or 
stays fixed, or is removed. This probabilistic interpretation of k 2 can be also extended 
to continuous random walk methods once we know the survival probability distribution 
function of a random walker in continuous space. 

In this letter, we obtain this survival probability distribution function of a random walker 
in continuous space by reinterpreting the weighting function in the previous modified WOS 



method [1 .1 for the LPBE (see the Appendix for more details). The survival probability of 



a random walker in a continuous and free diffusion region is given by JIT] : 

p(d) = dn/ smh(dK), (3) 

where d is the distance from the starting point in the diffusion region. Figure [I] shows 
this probability density function. We modify the WOS method to incorporate the survival 
probability to solve the LPBE via a continuous random walk method. This probability 
density combined with the WOS method is used to remove a random walker during the 



random walk by the acceptance-rejection method | 12]| . We generate a random number, r\ in 
[0, 1) when we perform a WOS step, and we compare 77 with p(d), the survival probability 
at d, the radius of WOS. If 77 > p(d), the random walker is removed at this WOS step. 
An estimate for the solution of the LPBE at xo, where random walkers start, is given by 

Sn- 
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SN = ^J2Mx n J. (4) 

iV i=l 

Here, N is the total number of random walkers, N s is the number of survived-and-absorbed 
random walkers, and X ni is the final position of the walker on the boundary when it is 
absorbed after rii WOS steps. 

In this method, like the WOS method, errors are due to both statistical sampling and 
the ^-absorption layer which captures random walkers near the boundary to terminate their 
random walk. However, the error from the 5-absorption layer can always be made smaller 
than the statistical error fIjT]]. For the same random walk, the estimate difference between 
using 5 and 5/10 gives a measure of the error due to the finite width of the ^-absorption layer. 
By adjusting 5 we can make the error from the absorption layer less than the statistical error. 
This means that if we increase the number of random walkers to decrease the statistical error, 
consequently we must reduce 5 and so increase the running time. 

In the following, we compare our simulation results with the analytic results for four 
problems, which were used as examples for the discrete random walk method 0. In all 
cases, the results are given as those normalized by the boundary condition ip Q , which is 
assumed sufficiently small for the LPBE to be valid. The number of random walks used for 
the solution at a point is 10 5 , and the absorption layer thickness is 5 = 10 -4 . The analytic 
results are shown with solid lines in Figs. 0, |] and |5| and our simulation results with 
circles. For the all four cases, our simulation results show excellent agreement. 

Our method has several features. First, it is easier to implement and will be faster 
than the other discretized methods, such as the discrete random walk method Q , the finite 
difference method [[HJ and the boundary-element method especially with complicated 



geometries. For a desired point, it takes only a few seconds to compute a solution with 
10 5 random walkers and 5 = 10~ 4 on a 550 Mhz PC. However, it is hard to compare 
to other methods because they compute solutions at all grid points. We can safely say 
that continuous Monte Carlo methods are more efficient when the solution is required only 
at relatively small number of points. Secondly, the accuracy and the running time of our 
method depends primarily on the number of statistical samples, and so it is naturally parallel. 



Thirdly, it is certain that our new method is faster than the old modified WOS method [i I 
because while some of our random walkers are removed during their random walk, in the old 
method all random walkers must complete their random walks to contribute to the solution 
according to their weightings. Also, in open boundary cases, like the three examples except 



the parallel plates, it is necessary to use a certain cut-off in the old modified WOS [|TT|] to 



kill random walkers, which will bias the results. As an example, in Table 1 in the case of 



the parallel plates, we compare our new method with the old modified WOS method, |TT 
We use the customary comparison method for Monte Carlo methods, the time consumption 
(or laboriousness) |15|]: t x D£, where t is the CPU time expended in calculating a single 



estimate and D£ is the variance of the estimates. The less laborious the algorithm, the more 
efficient it is. In Table 1, the time consumption (or laboriousness) of our algorithm is better 
than that of the old modified WOS method. Finally, our method is easy to extend to solve 



the LPBE with source terms ||11|| . That will be the subject of our upcoming research with 
biochemical applications. 

APPENDIX 
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In this appendix, we show how the weighting function in the old modified WOS 
method [[l l] can be interpreted as the survival probability distribution function. For sim- 
plicity, consider the LPBE in the old modified WOS method ]TT|. The solution at x in the 
domain can be expressed as follows 



1 



N 



(5) 



where 



Q° = 1, QT = Ql 



d'; 



n<— 1 



sinh(c? 



(6) 



Here, N is the total number of diffusing random walkers, i refers to ith random walker, X n . 
is the position where the ith random walker is absorbed in the (^-absorption layer after rij 
WOS steps, and d™* the radius of rijth WOS of the ith random walker. 

If we interpret Q™ 4 as a survival probability of ith random walker, X)i=i QT is the to- 
tal number of survived-and-absorbed random walkers. Furthermore, due to the property 
of probabilistic random sampling from the total random walkers, only the survived-and- 
absorbed random walkers can be regarded as contributors to the solution. This reinterpre- 
tation of the weighting function as the survival probability distribution function is a kind of 
the fractional sampling method, i.e. 'Russian Roullete', |12[ which has been used extensively 
in neutron transport and similar problems. 
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FIGURES 
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FIG. 1. The survival probability density function; d is the diffused distance of a random walker 
from the starting position, and k is the inverse Debye length. 
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r (in units of k l ) 

FIG. 2. The electric potential away from an charged infinite flat plate in an 0; the solid line 
is the analytic solution and the circles are the simulation results with I0 5 random walks and the 
absorption layer 5 = 10~ 4 . Here, r is the distance to the plate and n the inverse Debye length. 
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r (in units of k l ) 

FIG. 3. The electric potential in an electrolyte between two infinite charged parallel flat plates; 
the solid line is the analytic solution and the circles are the simulation results with 10 5 random 
walks and the absorption layer S = 10~ 4 . Here, r is the distance from the mid-point of the plates 
and k the inverse Debye length. 
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FIG. 4. The electric potential away from an infinitely long charged cylinder in an electrolyte; 
the solid line is the analytic solution and the circles are the simulation results with 10 5 random 
walks and the absorption layer 5 = 10 -4 . Here, r is the distance from the surface of the cylinder 
with unit radius and k the inverse Debye length. 
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r (in units of k" 1 ) 

FIG. 5. The electric potential away from the surface of a charged sphere in an electrolyte; the 
solid line is the analytic solution and the circles are the simulation results with 10 5 random walks 
and the absorption layer 5 = 10~ 4 . Here, r is the distance from the surface of the sphere with unit 
radius and n the inverse Debye length. 
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TABLES 



TABLE I. Time consumption comparison of our algorithm with the old modified WOS in the 
case of parallel plates at the mid-point; the variances are obtained from 100 independent runs, the 
number of random walks per run is 10 5 and the absorption layer, 5 = 10 -4 . 



method 


CPU time per run (sees) 


variance 


time consumption 


old method 
new method 


13.47 
2.97 


4.63 x 10" 7 
1.98 x 10" 6 


6.24 x 10" 6 
5.88 x 10~ 6 
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